Auxiliary Tool User Guide ************************************ .. note:: Want to quickly analyze results, plot data, or perform common data processing tasks after completing a DFT calculation with DSPAW? ``dspawpy`` (:bdg-info:`Python` :bdg-info-line:`>= 3.9`) is such a tool. It can be called programmatically (see example scripts below) and also provides a command-line interactive program. After following the tutorial :ref:`installation instructions `, you can use the interactive program by typing ``dspawpy`` and pressing Enter in the command line: .. code-block:: text ... loading dspawpy cli ... This is the dspawpy command-line interactive tool. Enjoy! ( ) _| | ___ _ _ _ _ _ _ _ _ _ _ _ /'_` |/',__)( '_`\ /'_` )( ) ( ) ( )( '_`\ ( ) ( ) ( (_| |\__, \| (_) )( (_| || \_/ \_/ || (_) )| (_) | \__,_)(____/| ,__/'`\__,_) \___x___/ | ,__/ \__, | | | | | ( )_| | (_) (_) \___/ Version: Installation Path ====================================== | 1: Update | 2: structure conversion | 3: Volumetric data processing | 4: Band structure calculation | 5: Density of States (DOS) data processing | 6: Joint display of band structure and density of states (DOS) | 7: Optical properties data processing | 8: NEB (Nudged Elastic Band) transition state calculation data processing | 9: phonon calculation data processing | 10: aimd ab initio molecular dynamics data processing | 11: Polarization data processing | 12: ZPE zero-point energy data processing | 13: Thermal correction energy of TS | | q: Quit ====================================== --> Enter a number and press Enter to select a function: Highlights: - Autocompletion: Works by pressing the Tab key, helping to quickly and correctly enter the required program arguments. - Multithreaded lazy loading: Loads modules in the background while waiting for user input, significantly reducing waiting time; loads only necessary modules, minimizing memory usage. Note: - When using on a remote server, the startup time may be longer due to poor disk I/O performance, potentially taking up to half a minute in extreme cases (directly related to the server's current disk I/O performance). If this is unacceptable, please install and use dspawpy on your own computer. - After typing `dspawpy` and pressing Enter, Python will first load built-in modules. Once this is complete, the prompt "... loading dspawpy cli ..." will appear, indicating the second stage (loading third-party dependencies) has begun. - After the second stage is completed, a welcome screen will be displayed, indicating that dspwapy has finished the initial loading and has entered the third stage. Subsequently, it will dynamically load the corresponding dependency libraries based on the selected functional modules, thereby minimizing waiting time. .. _install instruction: Installation and Updates ============================================================================== 1. On the HZW machine, dspawpy has been pre-installed. Activate the virtual environment using the following command to start using it: .. code-block:: shell source /data/hzwtech/profile/dspawpy.env 2. On other machines, please install dspawpy yourself (choose one of the following two methods): - Using mamba or conda, you can install the package from https://conda-forge.org/download/. .. code-block:: shell mamba install dspawpy -c conda-forge #conda install dspawpy -c conda-forge - Or, use pip3 (some operating systems may not have the executable `pip3`, in which case try `pip`) .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` pip :animate: fade-in-slide-down :color: info - pip3 is the package manager that comes with python3. - Linux and Mac usually come with Python 3 and pip3. - On Windows, open the Microsoft Store, search for Python, and install it. .. figure:: ../assets/python_ms_store.png :height: 300px :align: center Then open cmd or powershell to use pip. .. code-block:: shell pip3 install dspawpy For information on how to configure pip and conda mirror addresses to speed up the installation process, please refer to https://mirrors.tuna.tsinghua.edu.cn/help/pypi/ and https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/. If the installation still fails, try the mamba/conda installation method above. .. warning:: On clusters, due to permission issues, the `pip` in the public path may not support global installation of Python libraries. You must add the `--user` option after `pip install` to install them in your home directory under `~/.local/lib/python3.x/site-packages/`, where 3.x represents the Python interpreter version, and x can be any integer between 9 and 13. Python will prioritize loading dspawpy from the home directory, even if the version in the public environment is newer! Therefore, if you have previously installed dspawpy with ``--user`` and have forgotten to manually update the old version in your home directory, even after sourcing the public environment, you will not be able to call the dspawpy in the public environment. Instead, the old version will still be used, leading to some bugs. Therefore, **considering that the HZW cluster automatically updates dspawpy weekly, it is recommended not to install it redundantly in your home directory; delete any existing installations.** On other clusters, ensure that you manually update dspawpy in your home directory in a timely manner. If you prefer not to delete and update the `dspawpy` in your home directory, you can use the `-s` option when running your Python scripts to prevent importing `dspawpy` from your home directory: ``python -s your-script.py``. Update dspawpy --------------------------------------------------- To update dspawpy if it was installed with mamba/conda, use the following command: .. code-block:: shell mamba update dspawpy #conda update dspawpy If dspawpy was installed via pip: .. code-block:: shell pip install dspawpy -U # -U for upgrading to the latest version .. note:: If pip uses a domestic mirror site, it may fail to upgrade smoothly because the mirror site has not yet synchronized the latest version of `dspawpy `_ |dspawpy_version|. Please use the following command to tell pip to download and install from the official PyPI site: .. code-block:: shell pip install dspawpy -i https://pypi.org/simple --user -U # -i specifies the download address, --user installs for the current user only, and -U installs the latest version If you encounter errors related to **dspawpy** during runtime, first verify that you have correctly imported the latest version of dspawpy and check the installation path: .. code-block:: shell $ python3 # or python >>> import dspawpy >>> dspawpy.__version__ # will output the version number >>> dspawpy.__file__ # will output the installation path .. _structure conversion: structure structure conversion ============================================================================== To read structure information, use the ``read`` function; to write structure information to a file, use the ``write`` function; for quick structure conversion, use the ``convert`` function: .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: read(), write(), convert() :animate: fade-in-slide-down :color: info .. automodule:: dspawpy.io.structure :members: read, write, convert :noindex: See the :download:`2conversion.py <../../UserScripts/2conversion.py>` script for conversion: .. literalinclude:: ../../UserScripts/2conversion.py :language: python :linenos: The rules for setting several key parameters of the convert function are shown in the table below: .. list-table:: **dspawpy Supported IO format** :header-rows: 1 * - infmt (input file format) - infile (Input file name fuzzy match) - outfmt (output file format) - outfile (output filename fuzzy match) - Description * - h5 - `*.h5` - X - X - HDF5 files saved after DS-PAW calculations are completed * - json - `*.json` - json - `*.json` - json files saved after DS-PAW calculations are completed * - pdb - `*.pdb` - pdb - `*.pdb` - Protein Data Bank * - as - `*.as` - as - `*.as` - DS-PAW structure file containing atomic coordinates and other information * - hzw - `*.hzw` - hzw - `*.hzw` - DeviceStudio's default structure file * - xyz - `*.xyz` - xyz - `*.xyz` - Supports only single conformation of molecular structure when reading, and extended-xyz type trajectory files including unit cell when writing * - X - X - dump - `*.dump` - LAMMPS dump-type trajectory files * - X - `*.cif*/*.mcif*` - cif/mcif - `*.cif*/*.mcif*` - Crystallographic Information File * - X - `*POSCAR*/*CONTCAR*/*.vasp/CHGCAR*/LOCPOT*/vasprun*.xml*` - poscar - `*POSCAR*` - VASP files * - X - `*.cssr*` - cssr - `*.cssr*` - Crystal Structure Standard Representation * - X - `*.yaml/*.yml` - yaml/yml - `*.yaml/*.yml` - YAML Ain't Markup Language * - X - `*.xsf*` - xsf - `*.xsf*` - eXtended Structural Format * - X - `*rndstr.in*/*lat.in*/*bestsqs*` - mcsqs - `*rndstr.in*/*lat.in*/*bestsqs*` - Monte Carlo Special Quasirandom Structure * - X - `inp*.xml/*.in*/inp_*` - fleur-inpgen - `*.in*` - FLEUR structure file, requires the additional installation of the pymatgen-io-fleur library * - X - `*.res` - res - `*.res` - ShelX res structure file * - X - `*.config*/*.pwmat*` - pwmat - `*.config/*.pwmat` - PWmat files * - X - X - prismatic - `*prismatic*` - A file format used for STEM simulations * - X - `CTRL*` - X - X - Stuttgart LMTO-ASA files .. note:: - In the table above, `*` represents any character, and `X` indicates unsupported formats. - h5, json, pdb, xyz, dump, and CONTCAR formats support trajectory information consisting of multiple structures (common in structure optimization, NEB, or AIMD tasks) - The `in(out)fmt` parameter has higher priority than filename wildcard matching; for example, specifying `in(out)fmt='h5'` allows any filename, even `a.json`. - When writing structural information in `json` format, only visualization of NEB chain tasks is supported. See :ref:`observing neb chain` for details. - Structure information from DS-PAW output files such as neb.h5, phonon.h5, phonon.json, neb.json, and wannier.json is currently not readable. .. _ 电荷密度数据处理: Volumetric Data Processing ============================================================================== volumetricData Visualization ----------------------------------------------------- - See also :download:`3vis_vol.py <../../UserScripts/3vis_vol.py>`: .. literalinclude:: ../../UserScripts/3vis_vol.py :language: python :linenos: Drag the converted file :guilabel:`DS-PAW_rho.cube` into the **VESTA** software to visualize it: .. figure:: ../assets/1-3.png :height: 300px :align: center Differential volumetric data visualization ----------------------------------------------------- - See :download:`3dvol.py <../../UserScripts/3dvol.py>`: .. literalinclude:: ../../UserScripts/3dvol.py :language: python :linenos: The above script supports the processing of charge density differences in multi-component systems. As an example using a binary system, it generates the charge density difference file `delta_rho.cube` from AB.h5, A.h5, and B.h5. This file can be directly opened using **VESTA**. Volumetric data plane average ----------------------------------------------------- - See also :download:`3planar_ave.py <../../UserScripts/3planar_ave.py>`: .. literalinclude:: ../../UserScripts/3planar_ave.py :language: python :linenos: Processing the electrostatic potential file obtained from Section 3.3 of :doc:`/appliedcases` yields the following vacuum direction potential function curve: .. figure:: ../assets/us/3pot_ave.png :height: 450px :align: center .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: write_VESTA(), write_delta_rho_vesta(), average_along_axis() :animate: fade-in-slide-down :color: info - The ``write_VESTA`` function handles the visualization of volumetric data: .. automodule:: dspawpy.io.write :members: write_VESTA :noindex: volumetricData refers to physical quantities that vary with spatial position, such as charge density rho, potential energy function potential, localized charge density elf, partial charge density pcharge, and solvent-bound charge density rhoBound. This data is stored in the volumetricData type in DS-PAW. - The ``write_delta_rho_vesta`` function is responsible for handling the visualization of differential volumetricData: .. automodule:: dspawpy.io.write :members: write_delta_rho_vesta :noindex: - The ``average_along_axis`` function handles averaging volumetric data along a specific axis: .. automodule:: dspawpy.plot :members: average_along_axis :noindex: .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's possible that the program you are using (e.g., MobaXterm) is incompatible with the QT libraries. You should either switch programs (e.g., VSCode or the system's built-in terminal command line) or add the following code, starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') .. _ 能带数据处理: band data processing ============================================================================== .. note:: 1. The script calls get_band_data() to read data, and setting efermi=XX during data reading can shift the energy zero point to the specified value; setting zero_to_efermi=True can shift the energy zero point to the Fermi level in the read file. 2. When plotting using `pymatgen`'s `BSPlotter.get_plot()` in the script, you can set `zero_to_efermi=True` to shift the energy zero to the Fermi level. Due to a critical update in `pymatgen` on August 17, 2023, which changed the return object of the plotting function from `plt` to `axes`, subsequent scripts may become incompatible. Therefore, a conditional statement has been added to handle this in the relevant parts of the user's script. 3. For two-step band calculations, obtain the accurate Fermi level from the first-step self-consistent calculation (from the self-consistent `system.json`). If this fails, users can modify the energy zero point when calling `get_band_data` to read data, using the `efermi` parameter. For example: `band_data=get_band_data('band.h5',efermi=-1.5)` 4. When plotting, the script calls `BSPlotter.get_plot` from pymatgen. When the system is determined to be non-metallic, setting `zero_to_efermi` will consider the VBM as the Fermi level energy, rather than the Fermi level from the data file. Therefore, when the system is non-metallic, setting `zero_to_efermi=True` during data reading and setting `zero_to_efermi=True` during plotting will result in different plots. Running the Python script listed in this section, the program will determine whether the system is metallic. If it is a non-metallic system, you will be prompted to choose whether to shift the Fermi level to the zero energy point; please follow the prompts. Conventional Band Treatment ----------------------------------------------------------- See :download:`4bandplot.py <../../UserScripts/4bandplot.py>`: .. literalinclude:: ../../UserScripts/4bandplot.py :language: python :linenos: .. note:: For band structure calculations, an accurate Fermi level is required, which is obtained from the self-consistent calculation (from system.json). If the acquisition fails, users can modify the efermi parameter in the get_band_data function. Executing the code will generate a band structure plot similar to the following: .. figure:: ../assets/us/4bandplot.png :height: 300px :align: center The band is projected onto each element separately, with the size of the data points representing the element's contribution to the orbital. ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- See :download:`4bandplot_elt.py <../../UserScripts/4bandplot_elt.py>`: .. literalinclude:: ../../UserScripts/4bandplot_elt.py :language: python :linenos: .. note:: 1. To plot projected band structure data, use the BSPlotterProjected module. 2. Use the `get_elt_projected_plots` function in the `BSPlotterProjected` module to plot band diagrams with orbital contributions for each element. Executing the code will generate band plots similar to the following: .. figure:: ../assets/us/4bandplot_elt.png :height: 300px :align: center .. warning:: If you execute the script above by connecting to a remote server via SSH, and you encounter QT-related error messages, it's possible that the program you're using (such as MobaXterm) is incompatible with the QT libraries. You can either switch to a different program (e.g., VSCode or the system's built-in terminal command line), or add the following code, starting on the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Band projections onto different elements' different orbitals ------------------------------------------------------------------------------------------------------- Refer to :download:`4bandplot_elt_orbit.py <../../UserScripts/4bandplot_elt_orbit.py>`: .. literalinclude:: ../../UserScripts/4bandplot_elt_orbit.py :language: python :linenos: .. note:: 1. Use the `get_projected_plots_dots` method in the `BSPlotterProjected` module, which allows users to customize the band structure plots by specifying elements and orbitals (L) to be plotted. 2. For example, get_projected_plots_dots({'Mo': ['d'], 'S': ['s']}) plots the d-orbitals of Mo and the s-orbitals of S. Executing the code will generate a band structure plot similar to the following: .. figure:: ../assets/us/4bandplot_elt_orbit.png :height: 500px :align: center .. warning:: If you execute the above script by connecting to a remote server via SSH, and QT-related error messages appear, it may be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Projecting band structure onto different atomic orbitals -------------------------------------------------------- See :download:`4bandplot_patom_porbit.py <../../UserScripts/4bandplot_patom_porbit.py>`: .. literalinclude:: ../../UserScripts/4bandplot_patom_porbit.py :language: python :linenos: .. note:: 1. The `get_projected_plots_dots_patom_pmorb` function in the BSPlotterProjected module offers greater flexibility, allowing users to customize the band diagrams for specific atoms and orbitals. 2. Use `dictpa` to specify the atom, and `dictio` to specify the orbitals of that atom. 3. To superimpose projected components of some atoms or orbitals, specify the `sum_atoms` or `sum_morbs` parameters according to the documentation of the `get_projected_plots_dots_patom_pmorb` function. .. warning:: 1. If only a single orbital is selected and the orbital name has more than one letter (e.g., px, dxy, dxz), the get_projected_plots_dots_patom_pmorb function will raise an error. See `here `_ for details. Executing the code will generate band diagrams similar to the following: .. figure:: ../assets/us/4band_patom_porbit.png :height: 700px :align: center .. warning:: If you execute the above script by connecting to a remote server via SSH, and you encounter QT-related error messages, it's possible that the program you're using (e.g., MobaXterm) is incompatible with the QT libraries. Either switch to another program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Band unfolding processing --------------------------------------------------------- See :download:`4bandunfolding.py <../../UserScripts/4bandunfolding.py>`: .. literalinclude:: ../../UserScripts/4bandunfolding.py :language: python :linenos: Executing the code yields a band diagram similar to the following: .. figure:: ../assets/us/4bandunfolding.png :height: 400px :align: center .. warning:: .. warning:: This feature currently does not support setting the Fermi level of non-metallic materials as the zero-energy point (the default is the valence band top as the zero-energy point). .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it might be due to incompatibility between the program used (such as MobaXterm, etc.) and the QT library. You can either switch programs (e.g., VSCode or the system's built-in terminal) or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') band-compare band structure comparison figure processing --------------------------------------------------------- Plotting regular band structure and Wannier band structure on the same figure. Refer to :download:`4bandcompare.py <../../UserScripts/4bandcompare.py>`: .. literalinclude:: ../../UserScripts/4bandcompare.py :language: python :linenos: Executing the code will generate band comparison curves similar to the following: .. figure:: ../assets/phase4/wannier_band.png :height: 350px :align: center .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: get_band_data() :animate: fade-in-slide-down :color: info - The ``get_band_data`` function is responsible for reading band structure data as follows: .. automodule:: dspawpy.io.read :members: get_band_data :noindex: .. warning:: If you are running the above script by connecting to a remote server via SSH and encounter QT-related error messages, it may be due to incompatibility between the program you are using (such as MobaXterm, etc.) and the QT libraries. You can either switch to another program (such as VSCode or the system's built-in terminal), or add the following code, starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') .. _ 态密度数据处理: DOS Data Processing ============================================================================== Total Density of States ------------------------------------------------------- See :download:`5dosplot_total.py <../../UserScripts/5dosplot_total.py>`: .. literalinclude:: ../../UserScripts/5dosplot_total.py :language: python :linenos: .. note:: 1. Use the `get_dos_data` function to convert the `dos.h5` file obtained from DS-PAW calculations into a format supported by pymatgen. 2. Use the DosPlotter module to obtain the data from the DS-PAW calculated dos.h5 file. 3. The DosPlotter function can pass parameters: the `stack` parameter indicates whether to fill the DOS plots, and `zero_at_efermi` indicates whether to set the Fermi energy to zero in the DOS plot. Here, `stack=False` and `zero_at_efermi=False` are set. 4. Use `add_dos` in the `DosPlotter` module to add the DOS data. 5. Use the `get_plot` function in the `DosPlotter` module to plot the DOS. Executing the code will generate a density of states plot similar to the following: .. figure:: ../assets/us/5dos_total.png :height: 300px :align: center .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's likely that the program you're using (such as MobaXterm, etc.) is incompatible with the QT libraries. You can either switch programs (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Project Density of States onto different orbitals -------------------------------------------------------- See :download:`5dosplot_spd.py <../../UserScripts/5dosplot_spd.py>`: .. literalinclude:: ../../UserScripts/5dosplot_spd.py :language: python :linenos: .. note:: Use the `add_dos_dict` function in the `DosPlotter` module to obtain the projected density of states (DOS) data, and then use `get_spd_dos` to project the information onto spd orbitals. The code execution will produce a density of states plot similar to the following: .. figure:: ../assets/us/5dos_spd.png :height: 300px :align: center .. warning:: If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it might be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Projecting the density of states onto different elements --------------------------------------------------------- See also :download:`5dosplot_elt.py <../../UserScripts/5dosplot_elt.py>`: .. literalinclude:: ../../UserScripts/5dosplot_elt.py :language: python :linenos: .. note:: Use the `add_dos_dict` function in the DosPlotter module to obtain projected density of states data, then use `get_element_dos` to output the projected information according to different elements. The code execution will produce a density of states plot similar to the following: .. figure:: ../assets/us/5dos_elt.png :height: 300px :align: center .. warning:: If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it might be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Projecting the density of states onto different orbitals of different atoms ----------------------------------------------------------------------------------------------------------------- See :download:`5dosplot_atom_orbit.py <../../UserScripts/5dosplot_atom_orbit.py>`: .. literalinclude:: ../../UserScripts/5dosplot_atom_orbit.py :language: python :linenos: .. note:: 1. Use the `get_site_orbital_dos` function to extract the contribution of a specific atom and specific orbital from the DOS data. `dos_data.structure[0], Orbital(4)` represents obtaining the density of states for the dxy orbital of the first atom; the index in the `get_site_orbital_dos` function starts from 0. 2. Running this script and selecting the element and orbital as prompted will generate the corresponding density of states (DOS) plot. Executing the code will produce a density of states plot similar to the following: .. figure:: ../assets/us/5dos_atom_orbit.png :height: 300px :align: center .. warning:: If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it's likely due to incompatibility between the program you're using (e.g., MobaXterm) and the QT library. Either switch to a different program (such as VSCode or the system's built-in terminal command line), or add the following code to your Python script starting from the second line: .. code-block:: python import matplotlib matplotlib.use('agg') Projecting the density of states onto the split d-orbitals (t2g, eg) of different atoms --------------------------------------------------------------------------------------------------------------- See also :download:`5dosplot_t2g_eg.py <../../UserScripts/5dosplot_t2g_eg.py>`: .. literalinclude:: ../../UserScripts/5dosplot_t2g_eg.py :language: python :linenos: .. note:: 1. Use the `get_site_t2g_eg_resolved_dos` function to extract the t2g and eg orbital contributions for a specific atom from the DOS data. This retrieves the t2g and eg orbital contributions for the second atom. 2. Running this script and selecting an atom number as prompted will generate the corresponding density of states plot. Executing the code will generate a density of states plot similar to the following: .. figure:: ../assets/us/5dos_t2g_eg.png :height: 300px :align: center .. note:: If the element does not contain d orbitals, a blank image will be drawn. .. warning:: If you execute the script above by connecting to a remote server via SSH and encounter QT-related error messages, it's likely that the program you are using (such as MobaXterm) is incompatible with the QT library. You can either switch programs (e.g., VSCode or the system's built-in terminal command line) or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') d-centered analysis --------------------------------------------------------- Taking the Pb-slab system as an example, a d-band center analysis is performed on Pt atoms: See :download:`5center_dband.py <../../UserScripts/5center_dband.py>`: .. literalinclude:: ../../UserScripts/5center_dband.py :language: python :linenos: Executing the code yields results similar to the following: .. code-block:: python spin=1 -1.785319344084034 .. note:: Currently, only the d-orbital center averaged over all atoms is supported. Element-resolved, atom-projected, or other orbitals are not supported, nor is the selection of spin direction or energy range. The ``get_dos_data`` function is responsible for processing density of states data: .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: get_dos_data() :animate: fade-in-slide-down :color: info .. automodule:: dspawpy.io.read :members: get_dos_data :noindex: .. _ 能带和态密度共同显示: bandDos: Displaying Band Structure and Density of States Together ============================================================================== Using the Si system from the application tutorial as an example: Display band structure and density of states in a single figure. ------------------------------------------------------------------------------ See :download:`6bandDosplot.py <../../UserScripts/6bandDosplot.py>`: .. literalinclude:: ../../UserScripts/6bandDosplot.py :language: python :linenos: Executing the code yields a band density of states plot similar to the following: .. figure:: ../assets/us/6bandDos.png :height: 500px :align: center .. warning:: If you are connecting to a remote server via SSH and running the above script, and you encounter QT-related error messages, it's possible that the program you are using (such as MobaXterm) is incompatible with the QT libraries. You should either switch programs (e.g., VSCode or the system's built-in terminal command line) or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Display band structure and projected density of states on a single plot. ------------------------------------------------------------------------------ See :download:`6bandPdosplot.py <../../UserScripts/6bandPdosplot.py>`: .. literalinclude:: ../../UserScripts/6bandPdosplot.py :language: python :linenos: Executing the code yields a band-decomposed density of states plot similar to the following: .. figure:: ../assets/us/6bandPdos.png :height: 500px :align: center .. warning:: 1. Given projected band data, it will be projected along the element by default; given ordinary band data (or if the system contains more than 4 types of elements), it will not be projected and a warning will be output. 2. Given projected density of states (PDOS) data, projection along elements is also the default. You can switch to projection along orbitals, or no projection at all. For ordinary density of states (DOS) data and without disabling the DOS projection option BSDOSPlotter(dos_projection=None), the pymatgen plotting program will report an error, which is why a 6bandDosplot.py file was specifically prepared, as mentioned above. .. dspawpy还提供了一个优化后的band-compare能带对比图绘图函数 ``pltbd`` .. .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: pltbd() .. :animate: fade-in-slide-down .. :color: info .. .. automodule:: dspawpy.plot .. :members: pltbd .. :noindex: .. 新老函数绘图效果对比: .. .. figure:: ../assets/newbd.png .. :height: 500px .. :align: center .. .. figure:: ../assets/oldbd.png .. :height: 500px .. :align: center .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's likely due to incompatibility between the program you are using (e.g., MobaXterm) and the QT library. Either switch to a different program (such as VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') .. _ 光学性质数据处理: optical data processing ============================================================================== Using the :guilabel:`scf.h5` file obtained from a quick start calculation of the optical properties of the Si system as an example (Note: the output file name is the same as the task, task = scf; io.optical = true can calculate optical properties): Processing the reflectivity data, referring to :download:`7optical.py <../../UserScripts/7optical.py>`: .. literalinclude:: ../../UserScripts/7optical.py :language: python :linenos: .. note:: Reflectance is an optical property, and users can modify this keyword to "AbsorptionCoefficient", "ExtinctionCoefficient", or "RefractiveIndex" based on their needs, corresponding to the absorption coefficient, extinction coefficient, and refractive index, respectively. Executing the code will generate a curve showing the reflectance as a function of energy, similar to the following: .. figure:: ../assets/us/7optical.png :height: 300px :align: center .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: plot_optical() :animate: fade-in-slide-down :color: info .. automodule:: dspawpy.plot :members: plot_optical :noindex: .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it is likely that the program you are using (such as MobaXterm, etc.) is incompatible with the QT libraries. You can either switch to another program (such as VSCode or the system's built-in terminal command line) or add the following code to your Python script, starting from the second line: .. code-block:: python import matplotlib matplotlib.use('agg') neb data processing ============================================================================== Let's start with a quick introduction using the H diffusion on Pt(100) surface example: Generating intermediate configurations for input files ---------------------------------------------------------- - See :download:`8neb_interpolate_structures.py <../../UserScripts/8neb_interpolate_structures.py>`: .. literalinclude:: ../../UserScripts/8neb_interpolate_structures.py :language: python :linenos: .. note:: 1. Users can modify the number of interpolated points as needed. Setting it to 8 will generate a folder containing 8 structure files, with 6 intermediate configurations. 2. ``neb.linear_interpolate`` is a linear interpolation method. The ``pbc`` parameter, when set to ``True``, will lock the search for the shortest diffusion path. It defaults to ``False`` to increase user control, because 3. For example, if the initial fractional coordinate of an atom is 0.2 and the final state is 0.8. When pbc = True, the diffusion path will be forced to be 0.2 -> -0.2. When pbc = False, the user can make the program perform interpolation along the diffusion path 0.2 -> 0.8; if the shortest path is desired, manually change 0.8 to -0.2, thereby ensuring the program completes the initial guess of interpolation according to the user's intent. Plotting the energy barrier diagram ---------------------------------------------------------- neb.iniFin = true/false ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When neb.iniFin = true/false, you can use the path from the NEB calculation for barrier analysis (ensure that the initial and final state calculation files are in the NEB calculation path): - Refer to :download:`8neb_barrier_CubicSpline.py <../../UserScripts/8neb_barrier_CubicSpline_directory.py>`: .. literalinclude:: ../../UserScripts/8neb_barrier_CubicSpline_directory.py :language: python :linenos: .. note:: After running the above script, you can obtain a barrier curve similar to the following, with cubic spline interpolation: .. figure:: ./../assets/phase2/neb-barrier_cs.png :align: center :height: 500px For this specific example, the curve will exhibit an undesirable "dip" after cubic spline interpolation, which is inherent to the characteristics of the cubic spline interpolation algorithm. dspawpy internally calls `scipy's interpolation algorithms `_. Taking the cubic spline interpolation algorithm as an example in the script above, it is defined in the scipy documentation as: .. code-block:: python class scipy.interpolate.CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None) The keyword arguments include axis, bc_type, and extrapolate, whose specific meanings can be found in `scipy.interpolate.CubicSpline `_. We can specify the corresponding keyword arguments (axis, bc_type, extrapolate) in the plot_barrier function and pass them to the scipy.interpolate.CubicSpline class for processing. Here we use the script :download:`8neb_barrier.py <../../UserScripts/8neb_barrier.py>` to compare the curves plotted by interpolating with three algorithms: .. literalinclude:: ../../UserScripts/8neb_barrier.py :language: python :linenos: :emphasize-lines: 15, 21, 23 .. figure:: ./../assets/phase2/neb-barrier.png :height: 500px :align: center .. note:: 1. Choosing the appropriate interpolation algorithm is crucial for optimizing the final curve presentation. 2. In most cases, selecting the pchip (piecewise cubic Hermite interpolating polynomial) monotonic cubic spline interpolation algorithm will achieve good results, and it is also the default interpolation algorithm called. neb.iniFin = true ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When `neb.iniFin = true` is set, reading the `neb.h5/neb.json` files generated by the NEB calculation allows for a quick barrier analysis: - See :download:`8neb_barrier_CubicSpline.py <../../UserScripts/8neb_barrier_CubicSpline.py>`: .. literalinclude:: ../../UserScripts/8neb_barrier_CubicSpline.py :language: python :linenos: :emphasize-lines: 6 Processing the resulting barrier diagram is consistent with the previously read path. .. note:: 1. The energy stored in neb.h5 and neb.json files is TotalEnergy. If you need an accurate barrier value, it is recommended to process it by reading the NEB calculation path (taking TotalEnergy0). .. warning:: If you are connecting to a remote server via SSH to execute the script above and encounter QT-related error messages, it's possible that the program you are using (e.g., MobaXterm) is incompatible with the QT libraries. Either switch to a different program (such as VSCode or the system's built-in terminal command line), or add the following code starting on the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Processing Data for Transition State Calculations ---------------------------------------------------------- After NEB calculations, it is generally necessary to plot the energy barrier diagram and check the forces on each interpolated structure to ensure they are below a specified threshold. If the results are abnormal, the force and energy changes of each interpolated structure during the structure optimization process should also be checked to determine if they have truly "converged." These operations require at least three cycles. To simplify the process, we provide an all-in-one summary function ``summary``: - Refer to :download:`8neb_check_results.py <../../UserScripts/8neb_check_results.py>`: .. literalinclude:: ../../UserScripts/8neb_check_results.py :language: python :linenos: .. note:: 1. This script will print the energies and forces of each structure in a table, plot the energy barrier, and also plot the convergence of energy and forces for intermediate structures. 2. If ``neb.iniFin = false``, the user must copy the results file of the self-consistent calculation, either :guilabel:`scf.h5` or :guilabel:`system.json`, to the corresponding initial and final state subfolders. Otherwise, the program cannot read the energy and force information of the initial and final states and will exit with an error. 3. By default, the energy barrier plot is stored in the parent directory of the NEB calculation, and the energy and force convergence plots for each intermediate structure are stored in the respective subfolders. Executing the code will generate a table similar to the following, displaying the energy and force information for each NEB configuration: +-----+-------------+-----------------------+---------------------+---------------------------+ |Image|Force (eV/Å) |Reaction coordinate (Å)| Energy (eV) | Delta energy (eV) | +-----+-------------+-----------------------+---------------------+---------------------------+ |00 | 0.1803 | 0.0000 | -39637.0984 | 0.0000 | +-----+-------------+-----------------------+---------------------+---------------------------+ |01 | 0.0263 | 0.5428 | -39637.0186 | 0.0798 | +-----+-------------+-----------------------+---------------------+---------------------------+ |02 | 0.0248 | 1.0868 | -39636.8801 | 0.2183 | +-----+-------------+-----------------------+---------------------+---------------------------+ |03 | 0.2344 | 1.5884 | -39636.9984 | 0.1000 | +-----+-------------+-----------------------+---------------------+---------------------------+ |04 | 0.0141 | 2.0892 | -39637.0900 | 0.0084 | +-----+-------------+-----------------------+---------------------+---------------------------+ In addition to the energy barrier diagram, you can also obtain the energy and force convergence curves for each intermediate configuration (taking configuration 02 as an example). .. figure:: ./../assets/phase2/neb-energy.png :height: 500px :align: center .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it might be due to incompatibility between the program you are using (such as MobaXterm) and the QT library. Either switch to another program (such as VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') .. _observing neb chain: Observing the NEB Chain ---------------------------------------------------------- Here, the NEB chain refers to the geometric relationship between the interpolated structures (structure00.as, structure01.as, ...), rather than the changes of a single structure during the optimization process. - NEB calculations are computationally expensive, and observing the NEB chain helps to judge the convergence speed of the NEB calculation. Furthermore, after generating intermediate structures via interpolation, previewing the NEB chain is often necessary. These needs can be met using the :download:`8neb_visualize.py <../../UserScripts/8neb_visualize.py>` script: .. literalinclude:: ../../UserScripts/8neb_visualize.py :language: python :linenos: .. note:: 1. After this script generates the neb_movie*.json files, you can view them by opening the json file via ``Device Studio`` --> ``Simulator`` --> ``DS-PAW`` --> ``Analysis Plot``. 2. The `directory` parameter specifies the main path of the NEB calculation; the complete folder after the NEB calculation is finished must be provided. 3. This script supports processing ongoing (i.e., incomplete) NEB calculation files, allowing users to monitor the trajectory in real time. 4. The xyz file can be opened and viewed using OVITO software: Open the visualization interface via ``Device Studio`` --> ``Simulator`` --> ``OVITO``, and then drag and drop the xyz file. 5. Structure information reading priority: latestStructureXX.as > h5 > json; When ignorels is set to True, it first attempts to read data from h5, and if it fails, it reads from json. Calculate the inter-configuration distance ---------------------------------------------------------- - Refer to this script: :download:`8calc_dist.py <../../UserScripts/8calc_dist.py>`: .. literalinclude:: ../../UserScripts/8calc_dist.py :language: python :linenos: .. _neb_continuation_instructions: Continued calculation with neb -------------------------------------------------------- - To restart a NEB calculation, refer to :download:`8neb_restart.py <../../UserScripts/8neb_restart.py>`: .. literalinclude:: ../../UserScripts/8neb_restart.py :language: python :linenos: See :ref:`neb_continuation_instructions` for details. Energy and maximum atomic force variation trend during NEB calculation ----------------------------------------------------------------------------------------------------- - To view plots showing the energy and maximum atomic force trends during the NEB calculation, refer to :download:`8neb_energy_force_curves.py <../../UserScripts/8neb_energy_force_curves.py>`: .. literalinclude:: ../../UserScripts/8neb_energy_force_curves.py :language: python :linenos: Generates energy and force change trend charts: .. figure:: ./../assets/Energies.png :height: 500px :align: center .. figure:: ./../assets/MaxForce.png :height: 500px :align: center .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: write_neb_structures(), plot_barrier(), summary(), get_distance(), write_movie_json(), write_xyz(), restart() :animate: fade-in-slide-down :color: info - The `write_neb_structures` function is responsible for generating intermediate configurations: .. automodule:: dspawpy.diffusion.neb :members: write_neb_structures :noindex: - The `plot_barrier` function is responsible for plotting the energy barrier diagram: .. automodule:: dspawpy.diffusion.nebtools :members: plot_barrier :noindex: - The ``summary`` function is responsible for summarizing the NEB calculation task's documentation: .. automodule:: dspawpy.diffusion.nebtools :members: summary :noindex: - The ``get_distance`` function calculates the distance between two configurations: .. automodule:: dspawpy.diffusion.nebtools :members: get_distance :noindex: - The functions ``write_movie_json`` and ``write_xyz`` can write intermediate configurations to JSON or XYZ files: .. automodule:: dspawpy.diffusion.nebtools :members: write_movie_json, write_xyz :noindex: - The ``restart`` function is responsible for restarting the NEB calculation: .. automodule:: dspawpy.diffusion.nebtools :members: restart :noindex: - The ``monitor_force_energy`` function is responsible for plotting the energy and force changes during the NEB calculation: .. automodule:: dspawpy.diffusion.nebtools :members: monitor_force_energy :noindex: Phonon Data Processing ============================================================================== Using the example of a phonon band structure and density of states calculation for MgO, using :guilabel:`phonon.h5`: If phonopy is not installed, running the following script will result in the message ``no module named 'phonopy'``, but this does not affect the program's normal operation. Phonon band data processing ----------------------------------------------------- - Refer to :download:`9phonon_bandplot.py <../../UserScripts/9phonon_bandplot.py>`: .. literalinclude:: ../../UserScripts/9phonon_bandplot.py :language: python :linenos: Executing the code yields a phonon band structure curve similar to the following: .. figure:: ../assets/phase3/phonon-nonac.png :height: 500px :align: center .. warning:: If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it's likely due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code to your Python script starting from the second line: .. code-block:: python import matplotlib matplotlib.use('agg') Phonon Density of States Data Processing ---------------------------------------------------- - Refer to :download:`9phonon_dosplot.py <../../UserScripts/9phonon_dosplot.py>`: .. literalinclude:: ../../UserScripts/9phonon_dosplot.py :language: python :linenos: Executing the code yields a phonon density of states curve similar to the following: .. figure:: ../assets/us/9phonon_dosplot.png :height: 350px :align: center .. warning:: If you execute the script above by SSH connection to a remote server and encounter QT-related error messages, it's possible that the program you are using (such as MobaXterm) is incompatible with the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Phonon Thermodynamic Data Processing ----------------------------------------------------- Refer to :download:`9phonon_thermal.py <../../UserScripts/9phonon_thermal.py>`: .. literalinclude:: ../../UserScripts/9phonon_thermal.py :language: python :linenos: Executing the code yields phonon thermodynamic curves similar to the following: .. figure:: ../assets/us/9phonon.png :height: 500px :align: center .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: get_phonon_band_data(), get_phonon_dos_data(), plot_phonon_thermal() :animate: fade-in-slide-down :color: info - The ``get_phonon_band_data`` function is responsible for reading phonon band data: .. automodule:: dspawpy.io.read :members: get_phonon_band_data :noindex: - The ``get_phonon_dos_data`` function is responsible for reading the phonon density of states: .. automodule:: dspawpy.io.read :members: get_phonon_dos_data :noindex: - The ``plot_phonon_thermal`` function is responsible for plotting phonon thermodynamic properties: .. automodule:: dspawpy.plot :members: plot_phonon_thermal :noindex: .. warning:: If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's likely that the program you are using (e.g., MobaXterm) is incompatible with the QT libraries. You can either switch programs (e.g., VSCode or the system's built-in terminal) or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') .. _ 分子动力学模拟数据处理: aimd molecular dynamics simulation data processing ============================================================================== For a quick start, take the molecular dynamics simulation data of the :math:`H_{2}O` molecular system, for example, the :guilabel:`aimd.h5` file: Convert the trajectory file format to .xyz or .dump. ------------------------------------------------------ Read data from the HDF5 file output by AIMD and generate trajectory files. The generated .xyz or .dump files can be dragged and dropped into OVITO for visualization. You can open the OVITO visualization interface through Device Studio --> Simulator --> OVITO, and then drag and drop the .xyz or .dump files into OVITO. See :download:`10write_aimd_traj.py <../../UserScripts/10write_aimd_traj.py>`: .. literalinclude:: ../../UserScripts/10write_aimd_traj.py :language: python :linenos: Executing the code will generate trajectory files in .xyz and .dump formats, which can be opened with OVITO. For more details on structure file conversion, refer to :ref:`Structure Conversion` .. note:: OVITO and dspawpy do not support saving systems with non-orthogonal unit cells as dump files. Changes in energy, temperature, etc. curves during the dynamics process. ------------------------------------------------------------------------------------------------------ - Refer to :download:`10check_aimd_conv.py <../../UserScripts/10check_aimd_conv.py>`: .. literalinclude:: ../../UserScripts/10check_aimd_conv.py :language: python :linenos: Executing the code will generate the following combined graph: .. figure:: ../assets/phase2/aimd-joined.png :height: 800px :align: center .. warning:: If you execute the above script by SSH connection to a remote server and encounter QT-related error messages, it's possible that the program you are using (such as MobaXterm) is incompatible with the QT libraries. You can either switch programs (for example, VSCode or the system's built-in terminal command line), or add the following code starting from the second line of the Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Mean Squared Displacement (MSD) Analysis ------------------------------------------------------ - See :download:`10aimd_msd.py <../../UserScripts/10aimd_msd.py>`: .. literalinclude:: ../../UserScripts/10aimd_msd.py :language: python :linenos: Executing the code will generate an image similar to the following: .. figure:: ../assets/us/10MSD.png :height: 500px :align: center .. warning:: If you execute the above script by SSH connection to a remote server and encounter QT-related error messages, it might be due to incompatibility between the program you're using (like MobaXterm, etc.) and the QT libraries. Either switch to a different program (such as VSCode or the system's built-in terminal), or add the following code starting from the second line of the Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Root Mean Square Deviation (RMSD) Analysis ------------------------------------------------------ - See :download:`10aimd_rmsd.py <../../UserScripts/10aimd_rmsd.py>`: .. literalinclude:: ../../UserScripts/10aimd_rmsd.py :language: python :linenos: Executing the code will generate an image similar to the following: .. figure:: ../assets/us/10RMSD.png :height: 500px :align: center .. warning:: If you execute the above script by connecting to a remote server via SSH, and QT-related error messages appear, it may be due to incompatibility between the program used (e.g., MobaXterm) and the QT libraries. Either change the program (such as VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Analysis of Atomic Pair Distribution Functions or Radial Distribution Functions (RDFs) -------------------------------------------------------------------------------------------------------- - See :download:`10aimd_rdf.py <../../UserScripts/10aimd_rdf.py>` : .. literalinclude:: ../../UserScripts/10aimd_rdf.py :language: python :linenos: Executing the code will generate an image similar to the following: .. figure:: ../assets/us/10RDF.png :height: 500px :align: center - The statistical calculations involved in this section are complex; please refer to the function API for more details. .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: plot_aimd(), get_lagtime_msd(), plot_msd(), get_rs_rdfs(), plot_rdf(), get_lagtime_rmsd(), plot_rmsd() :animate: fade-in-slide-down :color: info - The ``plot_aimd`` function can be used to help check the convergence of key physical quantities during AIMD calculations: .. automodule:: dspawpy.plot :members: plot_aimd :noindex: - The ``get_*`` and ``plot_*`` functions are responsible for reading key physical quantities from the AIMD calculation process: .. automodule:: dspawpy.analysis.aimdtools :members: get_lagtime_msd, plot_msd, get_rs_rdfs, plot_rdf, get_lagtime_rmsd, plot_rmsd :noindex: - For manual data extraction and processing, refer to: .. code-block:: python :linenos: from dspawpy.io.read import get_sinfo from dspawpy.io.structure import read aimd_h5_files = ['aimd1.h5','aimd2.h5','aimd3.h5'] # Extract and merge data sequentially from multiple completed aimd.h5 files # Read data from multiple aimd.h5 files at once and create a list of pymatgen Structures pymatgen_structures = read(datafile=aimd_h5_files) # Or extract the arrays for i, df in enumerate(aimd_h5_files): # Get data from each aimd.h5 file sequentially Nstep, elements, positions, lattices, D_mag_fix = get_sinfo(df) .. warning:: If you connect to a remote server via SSH and execute the above script, and you encounter QT-related error messages, it's possible that the program you're using (such as MobaXterm) is incompatible with the QT libraries. Either change the program (for example, VSCode or the system's built-in terminal command line), or add the following code, starting on the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') Ferroelectric Polarization Data Processing ============================================================================== For a quick start, take the series of :guilabel:`scf.h5` files obtained from ferroelectric calculations on the :math:`HfO_{2}` system as an example: - See :download:`11Ferri.py <../../UserScripts/11Ferri.py>`: .. literalinclude:: ../../UserScripts/11Ferri.py :language: python :linenos: Executing the code will generate the following combined figure: .. figure:: ../assets/us/11pol.png :height: 500px :align: center The ferroelectric polarization values for the head and tail configurations can be found below: .. code-block:: python :linenos: from dspawpy.plot import plot_polarization_figure python plot_polarization_figure(directory='.', annotation=True, annotation_style=1) # Displays the ferroelectric polarization values for the initial and final configurations. The code will generate the following combined plot: .. figure:: ../assets/phase3/Ferri-Pola-anno1.png :height: 500px :align: center Alternatively, a second annotation style can be used: .. code-block:: python :linenos: from dspawpy.plot import plot_polarization_figure plot_polarization_figure(directory='.', annotation=True, annotation_style=2) # Displays the ferroelectric polarization values for the initial and final configurations. The code will generate the following combined plot: .. figure:: ../assets/phase3/Ferri-Pola-anno2.png :height: 500px :align: center .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: plot_polarization_figure() :animate: fade-in-slide-down :color: info - The `plot_polarization_figure` function is responsible for plotting the ferroelectric polarization figure: .. automodule:: dspawpy.plot :members: plot_polarization_figure :noindex: .. warning:: If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it may be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script: .. code-block:: python import matplotlib matplotlib.use('agg') .. _ ZPE零点振动能数据处理: Zero-Point Vibrational Energy Data Processing ============================================================================== Taking the :guilabel:`frequency.txt` file obtained from a quick start :math:`CO` system frequency calculation as an example, the zero-point vibrational energy is calculated based on the following formula: .. math:: ZPE=\sum_{i=1}^{3 N} \frac{h \nu_i}{2} where :math:`\nu_i` are the normal mode frequencies, :math:`h` is Planck's constant (:math:`6.626\times10^{-34} J\cdot s`), and :math:`N` is the number of atoms. - See :download:`12getZPE.py <../../UserScripts/12getZPE.py>`: .. literalinclude:: ../../UserScripts/12getZPE.py :language: python :linenos: The code execution results will be saved to the :guilabel:`ZPE.dat` file, and the file content is as follows: .. code-block:: text Data read from D:\quickstart\2.13\frequency.txt Frequency (meV) 284.840038 --> Zero-point energy, ZPE (eV): 0.142420019 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: getZPE() :animate: fade-in-slide-down :color: info - The ``getZPE`` function is responsible for calculating the zero-point vibrational energy: .. automodule:: dspawpy.io.utils :members: getZPE :noindex: .. _ TS的热校正能: TS Hot Calibration Energy ============================================================================== Contribution of the entropy change of the adsorbate to the energy ------------------------------------------------------------------------------------------------------- Calculation is based on the following formula: .. math:: \Delta S_{a d s}\left(0 \rightarrow T, P^0\right)=S_{v i b}=\sum_{i=1}^{3 N}\left[\frac{N_{\mathrm{A}} h \nu_i}{T\left(e^{h \nu_i / k_{\mathrm{B}} T}-1\right)}-R \ln \left(1-e^{-h \nu_i / k_{\mathrm{B}} T}\right)\right] Here, :math:`\Delta S_{a d s}` represents the entropy change of the adsorbate, calculated based on the harmonic approximation. :math:`S_{v i b}` represents the vibrational entropy. :math:`\nu_i` is the normal mode frequency, :math:`N_A` is Avogadro's constant (:math:`6.022\times10^{23} mol^{-1}`), :math:`h` is Planck's constant (:math:`6.626\times10^{-34} J\cdot s`), :math:`k_B` is the Boltzmann constant (:math:`1.38\times10^{-23} J\cdot K^{-1}`), :math:`R` is the ideal gas constant (:math:`8.314 J\cdot mol^{-1}\cdot K^{-1}`), and :math:`T` is the system temperature in units of :math:`K`. - See :download:`13getTSads.py <../../UserScripts/13getTSads.py>` for reference: .. literalinclude:: ../../UserScripts/13getTSads.py :language: python :linenos: The execution result will be saved to the :guilabel:`TS.dat` file, with the following content: .. code-block:: text Data read from D:\quickstart\2.13\frequency.txt Frequency (THz) 68.873994 --> Entropy contribution, T*S (eV): 4.7566990201851275e-06 Ideal gas entropy contribution to energy ----------------------------------------------------- Calculations are based on the following formula: .. math:: \begin{aligned} S(T, P) & =S\left(T, P^{\circ}\right)-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \\ & =S_{\text {trans }}+S_{\text {rot }}+S_{\text {elec }}+S_{\text {vib }}-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \end{aligned} Where: .. math:: S_{\text {trans }}=k_{\mathrm{B}}\left\{\ln \left[\left(\frac{2 \pi M k_{\mathrm{B}} T}{h^2}\right)^{3 / 2} \frac{k_{\mathrm{B}} T}{P^{\circ}}\right]+\frac{5}{2}\right\} .. math:: S_{\mathrm{rot}}= \begin{cases}0 & \text {, monatomic } \\ k_{\mathrm{B}}\left[\ln \left(\frac{8 \pi^2 I k_{\mathrm{B}} T}{\sigma h^2}\right)+1\right] & \text {, linear } \\ k_{\mathrm{B}}\left\{\ln \left[\frac{\sqrt{\pi I_A I_{\mathrm{B}} I_{\mathrm{C}}}}{\sigma}\left(\frac{8 \pi^2 k_B T}{h^2}\right)^{3 / 2}\right]+\frac{3}{2}\right\} & \text {, nonlinear }\end{cases} .. math:: S_{\text {elec }}=k_{\mathrm{B}} \ln [2 \times(\text { total spin })+1] .. math:: S_{\text {vib }}=k_{\mathrm{B}} \sum_i^{\text {vib DOF }}\left[\frac{\epsilon_i}{k_{\mathrm{B}} T\left(e^{\epsilon_i / k_{\mathrm{B}} T}-1\right)}-\ln \left(1-e^{-\epsilon_i / k_{\mathrm{B}} T}\right)\right] where :math:`I_A` to :math:`I_C` are the three principal moments of inertia for a non-linear molecule, :math:`I` is the degenerate moment of inertia for a linear molecule, and :math:`\sigma` is the symmetry number of the molecule. Furthermore, monatomic refers to a monatomic molecule, linear refers to a linear molecule, and nonlinear refers to a non-linear molecule. total spin is the total spin quantum number. vib DOF represents vibrational degrees of freedom. - Refer to the :download:`13getTSgas.py <../../UserScripts/13getTSgas.py>` script for processing: .. literalinclude:: ../../UserScripts/13getTSgas.py :language: python :linenos: .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: getTSads(), getTSgas() :animate: fade-in-slide-down :color: info - The ``getTSads`` function is responsible for calculating the contribution of adsorbate entropy change to the energy: .. automodule:: dspawpy.io.utils :members: getTSads :noindex: - The ``getTSgas`` function is responsible for calculating the contribution of ideal gas entropy change to energy: .. automodule:: dspawpy.io.utils :members: getTSgas :noindex: Appendix ============================================================================== - Quickly download all scripts by clicking :download:`UserScripts.zip <../dspawpy_proj/UserScripts.zip>` - `dspawpy Changelog `_